
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%This m-file generates Tables 7, 8 and Figure 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


clear all
clc
load BayesBeta_normal_1.mat

Age=Demo(:,1)';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%basic logit estimates
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Beta=mvnrnd(ones(S*N,1)*clogitbeta, Hessian(clogitbeta, Data));
Beta=Beta';
Beta=reshape(Beta,6,N,[]);

% Calculate the elasticity and marginal cost using the original coefficient
elas=elasticity(Beta(:,Age<22,:),Scenario);
elas_stderr=elasticity_stderr(Beta(:,Age<22,:),Scenario);
elas_prctile5=[elas_stderr(:,1) elas_stderr(:,3) elas_stderr(:,5)];
elas_prctile95=[elas_stderr(:,2) elas_stderr(:,4) elas_stderr(:,6)];
mc=Scenario(1,1)*(1+1/elas(1,1));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%1-mixture estimates
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


K=1;
D=3;

ncomp1=Bayesbeta(:,1:K)';                        %K by S matrix
ind1=Bayesbeta(:,K+1:K+N)';                      %N by S matrix
pvec1=Bayesbeta(:,K+N+1:K+N+K)';                 %K by S matrix
pvecmean1=mean(pvec1(:,5000:S),2);
Beta_test=Bayesbeta(:,K+N+K+1:K+N+K+L*N);   
Beta1=reshape(Beta_test',L,N,S);                 %L by N by S matrix  
clear Beta_test
Betamean1=reshape(mean(reshape(Beta1(:,:,5000:S),L*N,[]),2),L,N);
beta_test=Bayesbeta(:,K+N+K+L*N+1:K+N+K+L*N+L*D);
beta1=reshape(beta_test',D,L,S);                      %D by L by S matrix
clear beta_test
beta1mean=reshape(mean(reshape(beta1(:,:,5000:S),D*L,[]),2),D,L);
beta1std=reshape(std(reshape(beta1(:,:,5000:S),D*L,[]),0,2),D,L);
b_test=Bayesbeta(:,K+N+K+L*N+L*D+1:K+N+K+L*N+L*D+K*L);
b1=reshape(b_test',L,K,S);                             %L by K by S matrix
clear b_test
b1mean=reshape(mean(reshape(b1(:,:,5000:S),L*K,[]),2),L,K);
W_test=Bayesbeta(:,K+N+K+L*N+L*D+K*L+1:K+N+K+L*N+L*D+K*L+L*L*K);
W1=reshape(W_test',L,L,K,S);                          %L by L by K by S matrix
clear W_test
W1mean=reshape(mean(reshape(W1(:,:,:,5000:S),L*L*K,[]),2),L,L,K);

Beta1=Beta1(:,:,5001:S);
clear Bayesbeta

% Calculate the elasticity and marginal cost using the original coefficient
elas1=elasticity(Beta1(:,Age<22,:),Scenario);
elas1_stderr=elasticity_stderr(Beta1(:,Age<22,:),Scenario);
elas1_prctile5=[elas1_stderr(:,1) elas1_stderr(:,3) elas1_stderr(:,5)];
elas1_prctile95=[elas1_stderr(:,2) elas1_stderr(:,4) elas1_stderr(:,6)];
mc1=Scenario(1,1)*(1+1/elas1(1,1));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%3-mixture estimates
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load BayesBeta_normal_3.mat

Age=Demo(:,1)';

K=3;
D=3;

ncomp3=Bayesbeta(:,1:K)';                        %K by S matrix
ind3=Bayesbeta(:,K+1:K+N)';                      %N by S matrix
pvec3=Bayesbeta(:,K+N+1:K+N+K)';                 %K by S matrix
pvec3mean=mean(pvec3(:,5000:S),2);
Beta_test=Bayesbeta(:,K+N+K+1:K+N+K+L*N);   
Beta3=reshape(Beta_test',L,N,S);                 %L by N by S matrix  
clear Beta_test
Beta3mean=reshape(mean(reshape(Beta3(:,:,5000:S),L*N,[]),2),L,N);
beta_test=Bayesbeta(:,K+N+K+L*N+1:K+N+K+L*N+L*D);
beta3=reshape(beta_test',D,L,S);                      %D by L by S matrix
clear beta_test
beta3mean=reshape(mean(reshape(beta3(:,:,5000:S),D*L,[]),2),D,L);
beta3std=reshape(std(reshape(beta3(:,:,5000:S),D*L,[]),0,2),D,L);
beta3p5=reshape(prctile(reshape(beta3(:,:,5000:S),D*L,[]),5,2),D,L);
beta3p95=reshape(prctile(reshape(beta3(:,:,5000:S),D*L,[]),95,2),D,L);
b_test=Bayesbeta(:,K+N+K+L*N+L*D+1:K+N+K+L*N+L*D+K*L);
b3=reshape(b_test',L,K,S);                             %L by K by S matrix
clear b_test
b3mean=reshape(mean(reshape(b3(:,:,5000:S),L*K,[]),2),L,K);
W_test=Bayesbeta(:,K+N+K+L*N+L*D+K*L+1:K+N+K+L*N+L*D+K*L+L*L*K);
W3=reshape(W_test',L,L,K,S);                          %L by L by K by S matrix
clear W_test
W3mean=reshape(mean(reshape(W3(:,:,:,5000:S),L*L*K,[]),2),L,L,K);

Beta3=Beta3(:,:,5001:S);
clear Bayesbeta

% Calculate the elasticity and marginal cost using the original coefficient
elas3=elasticity(Beta3(:,Age<22,:),Scenario);
elas3_stderr=elasticity_stderr(Beta3(:,Age<22,:),Scenario);
elas3_prctile5=[elas3_stderr(:,1) elas3_stderr(:,3) elas3_stderr(:,5)];
elas3_prctile95=[elas3_stderr(:,2) elas3_stderr(:,4) elas3_stderr(:,6)];
mc3=Scenario(1,1)*(1+1/elas3(1,1));




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Posterior Mean of the beta (Table 7)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Table 7 (estimate, 5% percentile, and 95% percentile)
Table_7_estimate = beta3mean

Table_7_prctile5 = beta3p5

Table_7_prctile95 = beta3p95


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Price Elasticities (Table 8)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Table 8 (estimate, 5% percentile, and 95% percentile)
Table_8_estimate = [reshape(elas',[],1) reshape(elas1',[],1) reshape(elas3',[],1)]

Table_8_prctile5 = [reshape(elas_prctile5',[],1) reshape(elas1_prctile5',[],1) reshape(elas3_prctile5',[],1)]

Table_8_prctile95 = [reshape(elas_prctile95',[],1) reshape(elas1_prctile95',[],1) reshape(elas3_prctile95',[],1)]


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Fitted Densities for Random Coefficients (Figure 2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
eps=.0000001;
grid=-8:.2:10;
sim=1;
for s=5000:10:S
    density1(sim,1:length(grid),1:L)=0;
    density3(sim,1:length(grid),1:L)=0;
   for l=1:L
       test=inf*ones(length(grid),L);
       test(:,l)=grid';
       test2=test+eps;
       density1(sim,:,l)= ((mvncdf(test2,b1(:,1,s)',W1(:,:,1,s))-mvncdf(test,b1(:,1,s)',W1(:,:,1,s)))/eps)';
       for k=1:K
           density3(sim,:,l)= density3(sim,:,l)+(pvec3(k,s)...
               *(mvncdf(test2,b3(:,k,s)',W3(:,:,k,s))-mvncdf(test,b3(:,k,s)',W3(:,:,k,s)))/eps)';
       end
   end
   sim=sim+1;
end

for l=1:L
   mdensity1(l,:)=mean(density1(:,:,l));
   hdensity1(l,:)=mdensity1(l,:)+2*std(density1(:,:,l));
   ldensity1(l,:)=mdensity1(l,:)-2*std(density1(:,:,l));
   mdensity3(l,:)=mean(density3(:,:,l));
   hdensity3(l,:)=mdensity3(l,:)+2*std(density3(:,:,l));
   ldensity3(l,:)=mdensity3(l,:)-2*std(density3(:,:,l));
end

mdensityclogit(1:L,1:length(grid))=0;

for l=1:L
    test=1;
    while clogitbeta(l)>grid(test)
        test=test+1;
    end
    mdensityclogit(l,test)=max(density1(l,:));
end

figure(1)
Title_name = {'Legal Office'; 'Street Pirated Office'; 'Internet Pirated Office'; 'Price'; 'Download Time'; 'Update'};
for i=1:L
%    figure(i)
   subplot(3,2,i), plot(grid,mdensity1(i,:)','--k',grid,mdensity3(i,:)',':k')
   vline(clogitbeta(i),'-k','homo coef')
   legend('1 comp', '3 comp','fontsize',12)
   title(Title_name(i),'fontsize',12)
end

